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Abstract 



In recent work |14) , we have presented a novel approach for improv- 
ing particle filters for multi-target tracking. The suggested approach 
was based on drift homotopy for stochastic differential equations. Drift 
homotopy was used to design a Markov Chain Monte Carlo step which 
is appended to the particle filter and aims to bring the particle filter 
samples closer to the observations. In the current work, we present 
an alternative way to append a Markov Chain Monte Carlo step to a 
particle filter to bring the particle filter samples closer to the obser- 
vations. Both current and previous approaches stem from the general 
formulation of the filtering problem. We have used the currently pro- 
posed approach on the problem of multi-target tracking for both linear 
and nonlinear observation models. The numerical results show that 
the suggested approach can improve significantly the performance of a 
particle filter. 



Introduction 

Multi-target tracking is a central and difficult problem arising in many sci- 
entific and engineering applications including radar and signal processing, 
air traffic control and GPS navigation |llj . The tracking problem consists 
of computing the best estimate of the targets' trajectories based on noisy 
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measurements (observations). Several strategies have been developed for 
addressing the multi-target tracking problem, see e.g. [U [5j [151 E EH [HI 

m Ea [20]. 

As in our recent work [T3], in this paper we also focus on particle filter 
techniques 0[J5]. The popularity of the particle filter method has increased 
due to its flexibility to handle cases where the dynamic and observation 
models are non- linear and/or non-Gaussian. The particle filter approach is 
an importance sampling method which approximates the target distribution 
by a discrete set of weighted samples (particles). The weights of the samples 
are updated when observations become available in order to incorporate 
information from the observations. 

Despite the particle filter's flexibility, it is often found in practice that 
most samples will have a negligible weight with respect to the observation, 
in other words their corresponding contribution to the target distribution 
will be negligible. Therefore, one may resample the weights to create more 
copies of the samples with significant weights [8]. However, even with the 
resampling step, the particle filter might still need a lot of samples in order 
to approximate accurately the target distribution. Typically, a few samples 
dominate the weight distribution, while the rest of the samples are in sta- 
tistically insignificant regions. Thus, some authors (see e.g. [TJ [21]) have 
suggested the use of an extra step, after the resampling step, which can help 
move more samples in statistically significant regions. 

The extra step for the particle filter is a problem of conditional path sam- 
pling for stochastic differential equations (SDEs). In [17J, a new approach 
to conditional path sampling based on drift homotopy was presented. In 
that paper, it was also shown how the algorithm can be used to perform 
the extra step of a particle filter. In [T4] . we applied the conditional path 
sampling algorithm from [17J to perform the extra step of a particle filter 
for the problem of multi-target tracking. The numerical results in |14] sug- 
gested that the approach can improve significantly the performance of a 
particle filter for multi-target tracking. In the current work, we show yet 
another way of how to perform the extra step of a particle filter. Both the 
current approach and the one in [14] stem from the general formulation of 
the filtering problem. The details of the currently proposed implementation 
of the extra step for a single target are given in Section 11.31 and for multiple 
targets in Section 11.41 The relative merits of the proposed approach in this 
paper and the one proposed in [14] are briefly discussed in Section [3l A 
more detailed comparison will be presented elsewhere. 

To address the target-observation association problem we have used a 
simple Metropolis Monte Carlo algorithm which first appeared in [T3]. This 
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algorithm effects a probabilistic search of the space of possible associations 
to find the best target-observation association. Of course, one can use more 
sophisticated association algorithms (see [15] and references therein) but the 
Monte Carlo algorithm performed very well in the numerical experiments. 

The paper is organized as follows. Sections 11.11 and 11.21 provide a brief 
presentation of particle filters for single and multiple targets (more details 
can be found in [U EJ EJ US]), which will serve to highlight the versatility 
and drawbacks of this popular filtering method. Sect ions 1 1 . 3 1 and 1 1 . 4 1 demon- 
strate how one can use an extra step to improve the performance of particle 
filters for single and multiple targets. Section [2] presents numerical results 
for multi-target tracking for the cases of linear and nonlinear observation 
models. Finally, Section [3] contains a discussion of the results as well as 
directions for future work. 

1 Particle filtering 

Particle filters are a special case of sequential importance sampling methods. 
In Sections 11.11 and 11.21 we discuss the generic particle filter for a single and 
multiple targets respectively. In Sections 11.31 and [L4l we discuss the addition 
of an extra step to the generic particle filter for the cases of a single and 
multiple targets respectively. 

1.1 Generic particle filter for a single target 

Suppose that we are given an SDE system and that we also have access 
to noisy observations Zt x , • • • , Zt k of the state of the system at specified 
instants Ti,.., , Tk- The observations are functions of the state of the sys- 
tem, say given by = G(XT k ,Ck), where = 1,...,K are mutually 
independent random variables. For simplicity let us assume that the distri- 
bution of the observations admits a density g(Xx k , Zx k ), i.e., p(ZT k \XT k ) oc 
g(X Tk ,Z Tk ). 

The filtering problem consists of computing estimates of the conditional 
expectation E[f(XT k )\{ZT j }j = i\, i.e., the conditional expectation of the 
state of the system given the (noisy) observations. Equivalently, we are 
looking to compute the conditional density of the state of the system given 
the observations p(-XrJ{^r.,-}j = i)- There are several ways to compute this 
conditional density and the associated conditional expectation but for prac- 
tical applications they are rather expensive. 

Particle filters fall in the category of importance sampling methods. Be- 
cause computing averages with respect to the conditional density involves 
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the sampling of the conditional density which can be difficult, importance 
sampling methods proceed by sampling a reference density q{XT k \{ZT j }j=i) 
which can be easily sampled and then compute the weighted sample mean 

k i * p(X% \{Z T } k j=l ) 

E[f(X Tk )\{Z Tj }, =1 ] ^LW^p^Tj 



n=l 



or the related estimate 



,. , _ ^ n=i n ^ T " ) q(X^\{Z T .}^ 1 ) 

n=l q{xi\{z T M =1 ) 



E [f{ X T k )\{ZTj}j=i\ ~ (yn , {7 p ; 



where iV has been replaced by the approximation 

Particle filtering is a recursive implementation of the importance sampling 
approach. It is based on the recursion 

p(X Tk \{Z Tj }* =1 ) <xg(X Tk ,Z Tk )p(X Tk \{Z Tj }*ll), (2) 

where p(X Tk \{Z Tj })z\) = j 'p(^rJ^r fc _ 1 )p(^r fc _J{2T j }JZ 1 1 )dX Tfc _ 1 . (3) 

If we set 

?(X r J{ZrjJ =1 ) =p(^tJ{2t j }Jz 1 1 ), 
then from ([2]) we get 

Lfe 



p{x Tk \{z T .y j=1 ) 

, Y i, 7 Tfc \ « 9{XT k ,Z Tk 
q{X Tk \{Z T] Y j=1 ) 

The approximation in expression (Tjrj becomes 



En-1 f( x r)9( x r > z n) 
E[f(X T J\{Z Tj }U - ^ ^ I fc \ W 



From we see that if we can construct samples from the predictive distri- 



bution p(Xi< k \{Zi<- } „•_ J) then we can define the (normalized) weights Wp 
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g(X% ,Z T ) 

^tn fn= — -z — r, use them to weigh the samples and the weighted samples 

will be distributed according to the posterior distribution p(XT k \{Zr.}j = i )• 
In many applications, most samples will have a negligible weight with 
respect to the observation, so carrying them along does not contribute sig- 
nificantly to the conditional expectation estimate (this is the problem of 
degeneracy [9]). To create larger diversity one can resample the weights to 
create more copies of the samples with significant weights. The particle filter 
with resampling is summarized in the following algorithm due to Gordon et 
al. [8\. 

Particle filter for a single target 

1. Begin with N unweighted samples Xj, j from p(Xx k _ 1 \{ZTj }j=i)- 

2. Prediction: Generate N samples Xlf" from p{Xx k \XT k 

3. Update: Evaluate the weights 

(xtp k ,z Tk ) 



Wt = if 

k En=i9(X!p k ,Z Tk ) 

4. Resampling: Generate N independent uniform random variables 
{6 n }% =1 in (0, 1). For n = 1, . . . , N let X^ = X!} where 

zZ w k< e3 <J2 w k 

i=i i=i 

where j can range from 1 to N. 

5. Set k = k + 1 and proceed to Step 1. 

The particle filter algorithm is easy to implement and adapt for dif- 
ferent problems since the only part of the algorithm that depends on the 
specific dynamics of the problem is the prediction step. This has led to the 
particle filter algorithm's increased popularity [5J. However, even with the 
resampling step, the particle filter can still need a lot of samples in order to 
describe accurately the conditional density p{X r p h \ {Ztj }j=i ) • Snyder et al. 
[16j have shown how the particle filter can fail in simple high dimensional 
problems because one sample dominates the weight distribution. The rest 
of the samples are not in statistically significant regions. Even worse, as we 
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will show in the numerical results section, there are simple examples where 
not even one sample is in a statistically significant region. In the next sub- 
section we present how an extra step can be used to push samples closer to 
statistically significant regions. 

1.2 Generic particle filter for multiple targets 

Suppose that we have A = 1, . . . , A targets. Also, for notational simplicity, 
assume that the Ath target comes from the Ath observation. Even when this 
is not the case, we can relabel the observations to satisfy this assumption. 
The targets are assumed to evolve independently so that the observation 
weight of a sample of the vector of targets is the product of the individual 
observation weights of the targets [15]. The same is true for the transition 
density of the vector of targets between observations. We denote the vector 
of targets at observation Tjt by 

x n = (Xi,T k ,---,x A:Tk ) 

and the observation vector at Tk by 

Zr k = (Zi,T k ,- ■ -,Z A:Tk ). 

Also, we can have different observation weight densities g\, A = 1, . . . , A for 
different targets. However, in the numerical examples we have chosen the 
same observation weight density for all targets. 

Following |15| we can write the particle filter for the case of multiple 
targets as 

Particle filter for multiple targets 

1. Begin with A?" unweighted samples Xl^ f rom p( x T k _ 1 \{Zt 3 }j=i ) = 

2. Prediction: Generate N samples X'^ from 

A 
A=l 

3. Update: Evaluate the weights 

w n _ Yl\=l9\{X'x,T k i Z \T k ) 

En=l I! A=l 9\ ( X 'x,T k > Z \T k ) 
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4. Resampling: Generate N independent uniform random variables 
{6 n }% =1 in (0, 1). For n = 1, . . . , N let X" fc = X^ where 

E w k< 0J <E w k 

i=i i=i 
where j can range from 1 to N. 

5. Set k = k + 1 and proceed to Step 1. 

1.3 Particle filter with MCMC step for a single target 

Several authors (see e.g. [21]) have suggested the use of a MCMC step 
after the resampling step (Step 4) in order to move samples away from 
statistically insignificant regions. There are many possible ways to append 
an MCMC step after the resampling step in order to achieve that objective. 
The important point is that the MCMC step must preserve the conditional 
density p{Xj- k \ {ZTj}j = \)- in the current section we show that the MCMC 
step constitutes a case of conditional path sampling. 

We begin by noting that one can use the resampling step (Step 4) in the 
particle filter algorithm to create more copies not only of the good samples 
according to the observation, but also of the values (initial conditions) of the 
samples at the previous observation. These values are the ones who have 
evolved into good samples for the current observation (see more details in 
|21j). The motivation behind producing more copies of the pairs of initial 
and final conditions is to use the good initial conditions as starting points 
to produce statistically more significant samples according to the current 
observation. This process can be accomplished in two steps. First, Step 4 
of the particle filter algorithm is replaced by 

Resampling: Generate iV independent uniform random variables {9 n }^ = 
in (0, 1). For n = 1, ... ,N let (X^^X^) = (X^X^) where 

i=i i=i 

Also, through Bayes' rule [21] one can show that the posterior density 
p(XT k \{ZTj}j=i) is preserved if one samples from the density 

g(X Tk ,Z Tk )p(X Tk \X Tk ^) 
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where Xx k _ x are given by the modified resampling step. This is a problem 
of conditional path sampling for (continuous-time or discrete) stochastic 
systems. The important issue is to perform the necessary sampling efficiently 

mm- 

We are now in a position to present the particle filter with MCMC step 
algorithm 

Particle filter with MCMC step for a single target 

1. Begin with N unweighted samples Xj, from p(XT k _ x \{ZTj}jZi)- 

2. Prediction: Generate N samples Xlp from p{Xx k \XT k 

3. Update: Evaluate the weights 

w n _ yv r fc » lk, 



T — tTF 

k En=i9(X!p k ,Z Tk ) 

4. Resampling: Generate N independent uniform random variables 
{0"}£L ! in (0,1). For n = 1,...,N let (X^X^) = (X'^X^) 
where ^ 

zZ w k< ej <ll w k 

1=1 1=1 
where j can range from 1 to N. 

5. MCMC step: For n = 1,...,N choose a modified drift (possibly 
different for each n) . Construct a Markov chain for Yj! with stationary 
distribution 

g{Y\Z Tk ) V {Y n \X^_ x ) 

6. Set XS =YS . 

7. Set k = k + 1 and proceed to Step 1. 

Since the samples Xj, = Y^^ are constructed by starting from different 
sample paths, they are independent. Also, note that the samples Xj, are 
unweighted. However, we can still measure how well these samples approx- 
imate the posterior density by comparing the effective sample sizes of the 
particle filter with and without the MCMC step. For a collection of N 
samples the effective sample size essiTk) is defined by 

ess(T k ) = 
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where 



±J2(9(X? k ,Z Tk )-W k r and W k = ^ g(X^, Z Tk ). 

n=l n=l 

The effective sample size can be interpreted as that the N weighted samples 
are worth of ess(T k ) = j^-rz i.i.d. samples drawn from the target density, 

which in our case is the posterior density. By definition, ess{T k ) < N. If the 
samples have uniform weights, then ess(T k ) = N. On the other hand, if all 
samples but one have zero weights, then ess(T k ) = 1. 

1.4 Particle filter with MCMC step for multiple targets 

We discuss now the case of multiple, say A, targets. Instead of the observa- 
tions for a single target now we have a collection of observations for all the 
targets {Z r .}J =1 = {(Z*., . . . , )}* =1 . 

The particle filter with MCMC step for the case of multiple targets is 

Particle filter with MCMC step for multiple targets 

1. Begin with A" unweighted samples Xj, f rom P(XT k _ 1 \{ZT j }jZi) = 

2. Prediction: Generate A" samples Xlj^ from 

A 

p(X Tk \X Tk _ 1 ) = X[p{X KTk \X KTk _ 1 ). 

A=l 

3. Update: Evaluate the weights 

W n = U\=i9x( x '\,T^ Z \T k ) 

J2n=lU\=l9x(X'l Tk ,Zx,T k ) 

4. Resampling: Generate A" independent uniform random variables 
{<P}£U in (0,1). For n = 1,...,N let (X^X^) = {X'^X'^) 
where 

j'-i i 

zZ w k< ej <Y, w k 
i=i i=i 

where j can range from 1 to N. 
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5. MCMC step: For n = 1, . . . , N and A = 1, . . . , A choose a modified 
drift (possibly different for each n and each A). Construct a Markov 
chain for Yp with stationary distribution 

\[g x {Y^Z KTk )p x (YZ\Xl Tk J. 

A=l 

6. Set X™ = YS . 

7. Set k = k + 1 and proceed to Step 1. 

For a collection of N samples the effective sample size ess\(T k ) for A 
targets is 

ess\(T k ) 



N 



1 + cl k 



where 



W A , k \ 



N A 



n=l A=l 



^Y.(\\ax(xi Tk ,z XyTk )-w K , k f 

N A 

and W k , k = T7 E II 9A(Xl Tk ,Z x>Tk ). 



n=l A=l 



2 Numerical results 

We present numerical results for multi-target tracking using the particle 
filter with an MCMC step performed by hybrid Monte Carlo. We have syn- 
thesized tracks of targets moving on the xy plane using a 2D near constant 
velocity model [lj. At each time t we have a total of At targets and the 
evolution of the kth target (A = 1, . . . , At) is given by 



x-M = Ax Ajt _i + Bv Aii 

= [x\,t,x Xtt ,yx,t,yx,t] T , 



(5) 



where (x x< t, x\,t) and (y X j,y\,t) are the xy position and velocity of the fcth 
target at time t. The matrices A and B are given by 



1 T 

10 

1 T 

1 



and B 



T 2 /2 
T 








T 2 /2 
T 



(6) 
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where T is the time between observations. For the experiments we have set 
T = 1, i.e., noisy observations of the model are obtained at every step of 
the model © . The model noise V\,t is a collection of independent Gaussian 
random variables with covariance matrix H v defined as 







al 



(7) 



In the experiments we have a% = Uy = 1. Also, we have considered two 
possible cases for the observation model, one linear and one nonlinear. Due 
to the different possible combinations of targets to observations we use a 
different index m to denote the obsevations. Since we do not assume any 
clutter we have m = 1, . . . , A$. If the mth observation z mj j at time t comes 
from the A;th target we have 



y\,t 



+ w 



m,t 



(8) 



for the linear observation model and 



->m,t 



arctan 



y\,t 



+ w 



m,t 



(9) 



for the nonlinear observation model. As is usual in the literature, the non- 
linear observation model consists of the bearing 6 and range r of a target. 
The observation noise w m t is white and Gaussian with covariance matrix 



a 







obs,x 
a obs,y 



(10) 



for the linear observation model and 



a 2 e 
al 



(11) 



for the nonlinear observation model. For the numerical experiments with 



the linear observation model we chose a^ bs x 



a 



obs,y 



1. For the numerical 



experiments with the nonlinear observation model we chose a$ = 10 and 
a% = 1. These values make our example comparable in difficulty to examples 
appearing in the literature (see e.g. [151 EH 120])- 

The synthesized target tracks were created by specifying a certain sce- 
nario, to be detailed below, of surviving, newborn and disappearing targets. 
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According to this scenario we evolved the appropriate number of targets 
according to ([5]) and recorded the state of each target at each step. For 
the surviving targets we created an observation by using the state of the 
target in the observation model. Thus, for the linear observation model, 
the observations were created directly in xy space by perturbing the xy 
position of the target by (JS|) . For the nonlinear observation model, the ob- 
servations were created in bearing and range space 9,r by using Q. The 
perturbed bearing and range were transformed to xy space by the trans- 
formation x = rcos8, y = rsinO to create a position for the target in xy 
space. 

The newborn targets for the linear model were created in xy space di- 
rectly by sampling uniformly in [—100, 100]. Afterwards, the observations of 
the newborn targets were constructed by perturbing the x, y positions using 
([5]). The newborn targets for the nonlinear model were created in xy space 
by sampling uniformly in [—100, 100]. Afterwards, we transformed the x,y 
positions to the bearing and range space 9, r and perturbed the bearing and 
range according to ©. The perturbed bearing and range were again trans- 
formed back to xy space to create the position of the newborn target. Note 
that both observation models do not involve the velocities. The newborn 
target velocities were sampled uniformly in [—1, 1]. 

The number of targets at each observation instant is: Ao = 2, Ai = 2, 
A2 = 1, A3 = 2, A4 = 3, A5 = A6 = . . . = A200 = 4. So, for the majority 
of the steps we have 4 targets which makes the problem of tracking rather 
difficult. 

2.1 MCMC step 

For the nth sample, the density we have to sample for the linear model is 
A** 

II 9x(Zi jX)k , Z l,\,k)9y(zlx,ki Z 3,\,k)p( z \,k\ z \,k-l) 



A=l 



n( \ ( Z lXk z \Xk-l z 2,\,k~lT \ v x,\,k) 2 
exp v 1 ^ 

A=l v *• obs,x 

l^3,A,fc ^3,A ,fc-l ~ z i,X,k-l 1 ~ ~ V y,X,k) 



obs,y 



, Km) 2 . Kx, k ) 2 

2o-l + 2al 
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2.1.1 Hybrid Monte Carlo 

We chose to use Hybrid Monte Carlo to perform the sampling (any other 
MCMC method can be used) [9]. We present briefly the hybrid Monte Carlo 
(HMC) formulation that we have used to sample the conditional density for 
each target. 

Define the potential V{v™ k ,v™ k ) by 

^( u ",l,fej v y,l,k> ■■■ i v x,At k ,k> V y,At k ,k) = 

El \^l,\k z l,X,k-l Z 2Xk-\ L 2 V x,\,k> 
A=l obs,x 



+ 



{Z$Xk ~ z 3,X,k-l - z 2,X,k-l T - ^T V y,X,k) 2 



obs,y 



( v 2,X,k) 2 (- V y,X,kY i .-, 



Define the vectors t£ )fc = l . . . , < iAtfc ,J T and w£ fc = [u^ lifc , . . . , ^ Atfcjfc j 
where T is the transpose (not to be confuse with the interval between 
observations T used before). Consider v™ k ,Vy k as the position variables 
of a Hamiltonian system. We define the 2 At. -dimensional position vector 
q = [qi, q2\ T with q\ = v™ k and q2 = v™ k . To each of the position variables 
we associate a momentum variable and we write the Hamiltonian 

H(q,p)=V(q) + ^ 

where p = \pi,P2] T is the momentum vector. Thus, the momenta variables 
are Gaussian distributed random variables with mean zero and variance 1. 
The equations of motion for this Hamiltonian system are given by Hamilton's 
equations 

dqi dH dpi dH 

~r = and ~t~ = tor 1 = 1 ' 2 - 

dr opi dr oqi 

HMC proceeds by assigning initial conditions to the momenta variables 

T 

(through sampling from exp( — ^rp)), evolving the Hamiltonian system in 
fictitious time r for a given number of steps of size 5t and then using the 
solution of the system to perform a Metropolis accept/reject step (more de- 
tails in [9]). After the Metropolis step, the momenta values are discarded. 
The most popular method for solving the Hamiltonian system, which is the 
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Figure 1: Linear observation model. The solid lines denote the true target 
tracks, the crosses denote the observations and the dots the conditional 
expectation estimates from the improved particle filter. We have plotted the 
conditional expectation estimates every 5 observations to avoid cluttering 
in the figure. 

one we also used, is the Verlet leapfrog scheme. In our numerical imple- 
mentation, we did not attempt to optimize the performance of the HMC 
algorithm. For the sampling we used 100 Metropolis accept/reject steps 
and 1 HMC step of size St = 10 _1 to construct a trial path. 

For the nonlinear observation model we can use the same procedure 
as in the linear observation model to define a Hamiltonian system and its 
associated equations. We omit the details. 

2.2 Linear observations 

We start the presentation of our numerical experiments with results for the 
linear observation model (JH]) . Figures [T] and [2] show the evolution in the xy 
space of the true targets, the observations as well as the estimates of the 
improved particle filter. It is obvious from the figures that the improved 
particle filter follows accurately the targets and there is no ambiguity in the 
identification of the target tracks. 

The performance of the improved particle filter with 100 samples is com- 
pared to the performance of the generic particle filter with 120 samples in 
Figure [3] by monitoring the evolution in time of the RMS error per target. 
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Figure 2: Linear observation model. Detail of Figure [TJ 

The RMS error per target (RMSE) is denned with reference to the true 
target tracks by the formula 



RMSE(t) = 




where || ■ || is the norm of the position and velocity vector. Note that the state 
vector norm involves both positions and velocities even though the obser- 
vations use information only from the positions of a target, Xk,t is the true 
state vector for target k. E[x.k^\Zi, . . . , Zt] is the conditional expectation 
estimate calculated with the improved or generic particle filter depending 
on whose filter's performance we want to calculate. 

The improved particle filter has a computational overhead of the order of 
a few percent compared to the generic particle filter. We have thus used the 
generic particle filter with more samples than the improved particle filter. 
This additional number of samples more than accounts for the computational 
overhead of the improved particle filter. As can be seen in Figure [3] the 
generic particle filter's accuracy deteriorates quickly. On the other hand, 
the improved particle filter maintains an 0(1) RMS error per target for the 
entire tracking interval. The average value of the RMS error over the entire 
time interval of tracking is about 2.5 with standard deviation of about 0.5. 
For the generic particle filter, the average of the RMS error over the time 
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Figure 3: Linear observation model. Comparison of RMS error per target 
for the improved particle filter and the generic particle filter. 



interval of tracking is about 800 with standard deviation of about 760. 

Figure 0] compares the effective sample size for the generic particle filter 
and the improved particle filter. Because the number of samples is different 
for the two filters we have plotted the effective sample size as a percentage 
of the number of samples. We have to note that, after about 50 steps, 
the generic particle filter started producing observation weights (before the 
normalization) which were numerically zero. This makes the normalization 
impossible. In order to allow the generic particle filter to continue we chose 
at random one of the samples, since all of them are equally bad, and assigned 
all the weight to this sample. We did that for all the steps for which the 
observation weights were zero before the normalization. As a result, the 
effective sample size for the generic particle filter drops down to 1 sample 
after about 50 steps. Once the generic particle filter deviates from the true 
target tracks there is no mechanism to correct it. Also, we tried assigning 
equal weights to all the samples when the observation weights dropped to 
zero. This did not improve the generic particle's performance either. On 
the other hand, the improved particle filter maintains an effective sample 
size which is about 25% of the number of samples. 
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Figure 4: Linear observation model. Comparison of effective sample size for 
the improved particle filter and the generic particle filter. 
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Figure 5: Noninear observation model. The solid lines denote the true target 
tracks, the crosses denote the observations and the dots the conditional 
expectation estimates from the improved particle filter. We have plotted the 
conditional expectation estimates every 5 observations to avoid cluttering 
in the figure. 
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Figure 6: Nonlinear observation model. Detail of Figure [5j 



2.3 Nonlinear observations 

We continue with results for the nonlinear observation model ([9]) . Figures [5] 
and [6] show the evolution in the xy space of the true targets, the observations 
as well as the estimates of the improved particle filter. Again, as in the 
case of the linear observation model, the improved particle filter follows 
accurately the targets and there is no ambiguity in the identification of the 
target tracks. 

The case of the nonlinear observation model is much more difficult than 
the case of the linear observation model. The reason is that for the nonlinear 
observation model, the observation errors, though constant in bearing and 
range space, they become position dependent in xy space. In particular, 
when x and/or y are large, the observation errors can become rather large. 
This is easy to see by Taylor expanding the nonlinear transformation from 
bearing and range space to xy space around the true target values. Suppose 
that the true target bearing and range are #o> r o and its xy space position 
is xq = ro cos #0)2/0 = ros'm9o. Also, assume that the observation error in 
bearing and range space is, respectively, 59 and Sr. The xy position of a 
target that is perturbed by 89 and Sr in bearing and range space is (to first 
order) 

x = xq — yo59 — 5r cos 9q 
V = Uo + xq59 — 5r sin 9q- 
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Figure 7: Nonlinear observation model. Comparison of RMS error per target 
for the improved particle filter and the generic particle filter. 

Thus, the perturbation in xy space can be significant even if 59 and 5r are 
small. In our example we have ag = 10 -2 . So, when the true target x and 
y values become of the order of 10 3 as happens for some of the targets, the 
observation value in bearing and range space can be quite misleading as far 
as the xy space position of the target is concerned. As a result, even if one 
does a good job in following the observation in bearing and range space, the 
conditional expectation estimate of the xy space position can be inaccurate. 

With this in mind, we have used 200 samples for the improved particle 
filter and 220 samples for the generic particle filter. Again, the extra samples 
used for the generic particle filter more than account for the computational 
overhead of the improved particle filter. The performance of the improved 
particle filter is compared to the performance of the generic particle filter in 
Figure [7] by monitoring the evolution in time of the RMS error per target. 
The generic particle filter's accuracy again deteriorates rather quickly. The 
error for the improved particle filter is larger than in the linear observation 
model but never exceeds about 80 even after 200 steps when the targets have 
reached large values of x and/or y. The average value of the RMS error over 
the entire time interval of tracking is about 22 with standard deviation of 
about 21. For the generic particle filter, the average of the RMS error over 
the time interval of tracking is about 760 with standard deviation of about 
770. 
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Figure 8: Nonlinear observation model. Comparison of effective sample size 
for the improved particle filter and the generic particle filter. 

Figure [8] compares the effective sample size for the generic particle filter 
and the improved particle filter. After about 60 steps, the generic parti- 
cle filter, started producing observation weights (before the normalization) 
which were numerically zero. This makes the normalization impossible. In 
order to allow the generic particle filter to continue we chose at random 
one of the samples, since all of them are equally bad, and assigned all the 
weight to this sample. We did that for all the steps for which the observation 
weights were zero before the normalization. As a result, the effective sample 
size for the generic particle filter drops down to 1 sample after about 60 
steps. Once the generic particle filter deviates from the true target tracks 
there is no mechanism to correct it. Also, we tried assigning equal weights 
to all the samples when the observation weights dropped to zero. This did 
not improve the generic particle's performance either. On the other hand, 
the improved particle filter maintains an effective sample size which is about 
25% of the number of samples. 

3 Discussion 

We have presented an algorithm for multi-target tracking which builds on the 
existing particle filter methodology for multi-target tracking by appending 
an MCMC step after the particle filter resampling step. The purpose of the 
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addition of the MCMC step is to bring the samples closer to the observation. 
Even though the addition of an MCMC step for a particle filter has been 
proposed and used before [7], to the best of our knowledge, the currently 
proposed implementation of the MCMC step is novel (see also [21j for a 
related approach). 

We have tested the performance of the algorithm on the problem of 
tracking multiple targets evolving under the near constant velocity model [I] . 
We have examined two cases of observation models: i) a linear observations 
model involving the positions of the targets and ii) a nonlinear observation 
model involving the bearing and range of the targets. For both cases the 
proposed improved particle filter exhibited a significantly better performance 
than the generic particle filter. Since the improved particle filter requires 
more computations than the generic particle filter it is bound to be more 
expensive. However, the computational overhead of the improved particle 
filter is rather small, of the order of a few extra samples worth for the generic 
particle filter. 

In |14| we proposed another way of performing the extra MCMC step of 
a particle filter. That approach was based on modifying the drift of the dy- 
namic model and then accounting for the modification. In the current work, 
we use the original drift of the dynamic model without any modification. 
For the case of multi-target tracking with observations at every time step 
both algorithms perform equally well. Thus, at first sight it would appear 
that there is no need for the extra complication of modifying the drift of 
the dynamic model and then accounting for the modification as was done in 
|14j . However, in cases where there are only sparse observations, the sam- 
pling of the conditional density needed for the extra step can be much more 
difficult (and consequently expensive) for the original dynamic model than 
for the modified dynamic model. With this in mind, the approach in p3] 
seems to have wider applicability. A detailed comparison of the algorithm 
proposed in the current work and the one proposed in [14] will be presented 
in a future publication. 
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